function main

clearvars;

% =========================================================================
% 使用预生成好的图像作为介质进行超声和光声仿真
% =========================================================================

%% define the source and destination paths
source_folder_path = 'D:/project/dataset/dataset_with_noise/CARTILAGE_PLATES_Dataset';
dst_us_path = 'D:/project/dataset/ay_net_dataset/US_dataset';
dst_pa_path = 'D:/project/dataset/ay_net_dataset/PA_dataset';

%% define some constants to be used on the script
simulation_target_type = 'CARTILAGE';  % simulation_target的可选参数为 'BONE'，'VESSEL', 'CARTILAGE'
number_of_plane_waves = 1;          % number of plane waves in CPWC sequence  激发的平面波数量，先用1

%% define the ultrasound probe.
%% 2.5MHZ 探头参数==================%
%% PA&US shared reconstruction parameters =========%
params.Nelements = 64;
params.fs = 40e6;   % [Hz]
params.fc = 2.5e6;  % [Hz]
params.bandwidth = 90;    % [percent]
params.pitch = 0.300e-3;  % pitch == spacing [mm]
params.width = 0.254e-3;  %  [mm]
params.radius = inf;   % 线阵
params.c = 1500;  %[m/s]
params.fnumber = [];
%[X, Z] = meshgrid(x, z); % 生成网格点坐标[mm]

% %% 10MHZ 探头参数==================%
% %% PA&US shared reconstruction parameters =========%
% params.Nelements = 64;
% params.fs = 40e6;   % [Hz]
% params.fc = 10e6;  % [Hz]
% params.bandwidth = 90;    % [percent]
% params.pitch = 0.15e-3;  % pitch == spacing [mm]
% params.width = 0.127e-3;  %  [mm]
% params.radius = inf;   % 线阵
% params.c = 1500;  %[m/s]
% params.fnumber = [];
% xlen = params.pitch * 63 + params.width;
% x = (-xlen/2):0.025e-3:(xlen/2);
% z = 0:0.025e-3:22.5e-3;  % exp 60mm = 60e-3
% [X, Z] = meshgrid(x, z); % 生成网格点坐标[mm]


% 获取目录下的所有文件和子目录信息
dirInfo = dir(source_folder_path);
% 过滤出文件（排除 . 和 ..）
files = dirInfo([dirInfo.isdir] == 0);
if isempty(files)
    disp('No files found in the specified source_folder_path.');
end
% 提取文件名
fileNames = {files.name};


dst_dirInfo = dir(dst_pa_path);
% 过滤出文件（排除 . 和 ..）
dst_files = dst_dirInfo([dst_dirInfo.isdir] == 0);
dst_fileNames = {dst_files.name};


% 依次遍历 fileNames 中的每个文件名，拼接出文件完整路径并打印出来
    for k = 1:length(fileNames)
        if k < 629
            continue;
        end

        [~, baseName, ext] = fileparts(fileNames{k});

%         has_processed = false; 
%         for m = 1 : length(dst_fileNames)
%             [~, dst_baseName, ext] = fileparts(dst_fileNames{m});
%             if contains(dst_baseName, baseName)
%                 has_processed = true;
%                 break;
%             end
%         end
%         if has_processed
%             disp(fprintf('---- %s 已处理，跳过 ----', baseName));
%             disp(fprintf('---- k === %d ----', k));
%             continue;
%         end

        target_path = fullfile(source_folder_path, fileNames{k});
        [us_rf_data, us_img, pa_rf_data, pa_img] = us_pa_simulation(target_path, simulation_target_type, params, number_of_plane_waves);

        %%将探头收集到的超声RF数据存盘
        us_rf_data_path = fullfile(dst_us_path, [baseName '_us_rf.mat']);
        rf_data = us_rf_data.rf_data;
        acoustic_impedance = us_rf_data.acoustic_impedance;
        save(us_rf_data_path, 'rf_data', 'acoustic_impedance');
        % 将重建图像的数据范围缩放到 [0, 255] 以便于保存为图像
        scaled_us_image = mat2gray(us_img);
        % 使用 imwrite 函数保存图像
        us_image_path = fullfile(dst_us_path, [baseName '_us_das.png']);
        imwrite(uint8(scaled_us_image * 255), us_image_path);
        % 将源图像拷贝到target目录作为ground truth        
        us_image_gt_path = fullfile(dst_us_path, [baseName '_us_gt.png']);
        copyfile(target_path, us_image_gt_path);

        %%%  将探头收集到的光声RF数据存盘以备后续网络使用
        pa_rf_data_path = fullfile(dst_pa_path, [baseName '_pa_rf.mat']);
        rf_data = pa_rf_data.rf_data;
        optical_absorption = pa_rf_data.optical_absorption;
        save(pa_rf_data_path, 'rf_data', 'optical_absorption');
        % 将重建图像的数据范围缩放到 [0, 255] 以便于保存为图像
        scaled_pa_image = mat2gray(pa_img);
        % 使用 imwrite 函数保存图像
        pa_image_path = fullfile(dst_pa_path, [baseName '_pa_das.png']);
        imwrite(uint8(scaled_pa_image * 255), pa_image_path);
        pa_image_gt_path = fullfile(dst_pa_path, [baseName '_pa_gt.png']);
        copyfile(target_path, pa_image_gt_path);

        disp(fprintf('---- 成功保存%d组超声和光声数据到指定目录 ----', k));
    end
end 


%%% 对一个目标对象进行超声和光声仿真，目标对象的输入形式为图像，
%%% 在超声仿真时此图像表示为介质的密度和声速分布，光声仿真时此图像表示为初始声压分布（声源）
%%% 输出为仿真得到的超声RF data（us_rf_data），DAS重建得到的超声图像（us_img），光声RF data（pa_rf_data），DAS重建得到的光声图像（pa_img） 
%%  TODO：后续可能增加输入参数，如图像重建方法
%%
%% target_path: 图像文件（目标）的完整路径
%% target_type: 目标的类型，可选值为 'VESSEL' ,'BONE'或 'CARTILAGE'
%% params: 换能器参数，超声时用于发射和接收，光声时仅用于接收
%% num_plane_waves:  plane waves number in CPWC sequence   激发的平面波数量，先用1
function [us_rf_data, us_img, pa_rf_data, pa_img] = us_pa_simulation(target_path, target_type, params, num_plane_waves)
    % 打印文件完整路径
    disp(['simulation target:', target_path]);
    %% define some constants to be used on the script
    cycles = 2;       % number of cycles in pulse    激发信号含几个周期
    c0 = 1540;      % medium speed of sound [m/s]  超声在水中传播的速度
    rho0 = 1020;    % medium density [kg/m3]  水的密度
    F_number = 1.7; % F number for CPWC sequence (i.e. maximum angle)  含义是什么？
    PML_size = 20;         % size of the PML in grid points
    N = num_plane_waves;   % number of plane waves in CPWC sequence  激发的平面波数量，先用1


    xlen = params.pitch * (params.Nelements - 1) + params.width;
    x_imaging_axis = (-xlen/2):0.05e-3:(xlen/2);
    z_imaging_axis = 0:0.05e-3:40e-3;  % exp 60mm = 60e-3
    %[X, Z] = meshgrid(x, z); % 生成网格点坐标[mm]

    %% define the computational grid as a uff.linear_scan strcuture. 
    %% We set different resolution options depending on frequency reference speed of sound.
    f0 = params.fc;       % pulse center frequency [Hz]  激发信号的中心频率
    f_max = 1.2*f0;
    lambda_min = c0/f_max;

    % mesh resolution, choose one
    mesh_resolution = 'element4'; 
    switch mesh_resolution
        case 'element2' % around 50 sec per wave
            dx=params.pitch/2;                            % 2 elements per pitch 
        case 'element4' % around 6min sec per wave
            dx=params.pitch/4;                            % 4 elements per pitch 
        otherwise
            error('Not a valid option');
    end

    p0 = loadImage(target_path); 
    %%%p0 = imbinarize(p1);  %%% 仅用于复杂结构的血管VESSEL

    % mesh size
    Nx = size(p0, 2);
    Nz = size(p0, 1);
    grid_width = Nx*dx;    %% 真实的x方向视野, 大约2cm
    grid_depth = Nz*dx;    %% 真实的z方向视野，大约4cm
    % domain=uff.linear_scan('x_axis', linspace(-grid_width/2,grid_width/2,Nx).', 'z_axis', linspace(0,grid_depth,Nz).');

    kgrid = kWaveGrid(Nz, dx, Nx, dx); %% tricky: 与输入的图像矩阵大小须相同
    gridz = kgrid.x-kgrid.x_vec(1);  %% gridz和gridx仅用于显示坐标轴上的刻度值
    gridx = kgrid.y;

%     figure;
%     imagesc(gridx(1,:)*1e3, gridz(:,1)*1e3, p0); colormap gray; colorbar; axis equal tight;
%     xlabel('x [mm]');
%     ylabel('z [mm]');
%     title('original image');

    %% define the medium based by setting the sound speed and density in every pixel of the uff.scan.
    % transparent background
    medium.sound_speed = c0*ones(Nz, Nx);   % sound speed [m/s]
    medium.density = rho0.*ones(Nz, Nx);      % density [kg/m3]

% %     if strcmp(target_type, 'BONE') == 1
% %         disp('simulation target Bone');
% %         medium.sound_speed(p0<1) = 3360;       % 骨头的声速，sound speed [m/s]这儿的设置是否准确？ 颅骨声速3360，血液声速1570，人体软组织声速1500
% %         medium.density(p0<1) = 1658;           % 骨头的密度， density [kg/m3]  颅骨密度1658，血液密度1055，人体软组织1016
% %         acoustic_impedance = 5.57;
% %         optical_absorption = 0.8;                %% 680nm下，光声信号强度 血液>硬骨(HAP)>软骨（Cartilage）>水
% %     elseif strcmp(target_type, 'VESSEL') == 1
% %         disp('simulation target Vessel');
% %         medium.sound_speed(p0<1) = 1570;       % 血液的声速，sound speed [m/s]这儿的设置是否准确？ 颅骨声速3360，血液声速1570，人体软组织声速1500
% %         medium.density(p0<1) = 1055;           % 血液的密度， density [kg/m3]  颅骨密度1658，血液密度1055，人体软组织1016
% %         acoustic_impedance = 1.656;
% %         optical_absorption = 1;                %% 680nm下，光声信号强度 血液>硬骨(HAP)>软骨（Cartilage）>水
% %     elseif strcmp(target_type, 'CARTILAGE') == 1
% %         disp('simulation target Cartilage');
% %         medium.sound_speed(p0<1) = 1627;       % 软骨的声速，sound speed [m/s]这儿的设置是否准确？ 颅骨声速3360，血液声速1570，人体软组织声速1500
% %         medium.density(p0<1) = 1100;           % 软骨的密度， density [kg/m3]  颅骨密度1658，血液密度1055，人体软组织1016， 900， 1000， 1100
% %         acoustic_impedance = 1.87;
% %         optical_absorption = 0.4;               %% 680nm下，光声信号强度 血液>硬骨(HAP)>软骨（Cartilage）>水
% %     else
% %         disp('Invalid simulation_target value');
% %     end

    p0((p0 > 0.1) & (p0 ~= 1)) = 0;                             %%%图像插值后图像不是真正的二值图像，而是存在灰度区域，在仿真时将这些灰度区域也设置为样本区域，以保持与不添加扰动时的结果接近
    if strcmp(target_type, 'BONE') == 1
        disp('simulation target Bone');
        medium.sound_speed(p0<1) = 3360 * (1 - p0(p0<1));       % 骨头的声速，sound speed [m/s]这儿的设置是否准确？ 颅骨声速3360，血液声速1570，人体软组织声速1500
        medium.density(p0<1) = 1658  * (1 - p0(p0<1));           % 骨头的密度， density [kg/m3]  颅骨密度1658，血液密度1055，人体软组织1016
        acoustic_impedance = 5.57;
        optical_absorption = 0.8;                %% 680nm下，光声信号强度 血液>硬骨(HAP)>软骨（Cartilage）>水
    elseif strcmp(target_type, 'VESSEL') == 1
        disp('simulation target Vessel');
        medium.sound_speed(p0<1) = 1570  * (1 - p0(p0<1));       % 血液的声速，sound speed [m/s]这儿的设置是否准确？ 颅骨声速3360，血液声速1570，人体软组织声速1500
        medium.density(p0<1) = 1055  * (1 - p0(p0<1));           % 血液的密度， density [kg/m3]  颅骨密度1658，血液密度1055，人体软组织1016
        acoustic_impedance = 1.656;
        optical_absorption = 1;                %% 680nm下，光声信号强度 血液>硬骨(HAP)>软骨（Cartilage）>水
    elseif strcmp(target_type, 'CARTILAGE') == 1
        disp('simulation target Cartilage');
        medium.sound_speed(p0<1) = 1627  * (1 - p0(p0<1));       % 软骨的声速，sound speed [m/s]这儿的设置是否准确？ 颅骨声速3360，血液声速1570，人体软组织声速1500
        medium.density(p0<1) = 1100  * (1 - p0(p0<1));           % 软骨的密度， density [kg/m3]  颅骨密度1658，血液密度1055，人体软组织1016， 900， 1000， 1100
        acoustic_impedance = 1.87;
        optical_absorption = 0.4;               %% 680nm下，光声信号强度 血液>硬骨(HAP)>软骨（Cartilage）>水
    else
        disp('Invalid simulation_target value');
    end
    
    % attenuation
    medium.alpha_coeff = 0.3;  % [dB/(MHz^y cm)]
    medium.alpha_power = 1.5;
    
% % %     % show physical map: speed of sound and density
% % %     figure;
% % %     subplot(1,2,1);
% % %     imagesc(gridx(1,:)*1e3, gridz(:,1)*1e3, medium.sound_speed); colormap hot; colorbar; axis equal tight;
% % %     xlabel('x [mm]');
% % %     ylabel('z [mm]');
% % %     title('c_0 [m/s]');
% % %     subplot(1,2,2);
% % %     plot(medium.sound_speed(:,134));
% % %     
% % %     figure;
% % %     subplot(1,2,1);
% % %     imagesc(gridx(1,:)*1e3, gridz(:,1)*1e3, medium.density); colormap hot; colorbar; axis equal tight;
% % %     xlabel('x [mm]');
% % %     ylabel('z [mm]');
% % %     title('\rho [kg/m^3]');
% % %     subplot(1,2,2);
% % %     plot(medium.density(:,134));
    
    %% define the time vector depending on the CFL number, the size of the domain and the mean speed of sound.
    cfl = 0.3;  %% YCS：这个参数的含义
    t_end = 2*sqrt(grid_depth.^2+grid_depth.^2)/mean(medium.sound_speed(:));  %% 超声沿网格对角线往返的时间
    kgrid.makeTime(medium.sound_speed, cfl, t_end);
    kgrid.Nt=round(kgrid.Nt/1.5);%% 可以删除
    
    %% define a sequence of plane-waves，多个角度平面波的代码先保留
    % alpha_max=1/2/F_number;                       % maximum angle span [rad]
    alpha_max = 10/180*pi;
    if N>1
        angles=linspace(-alpha_max,alpha_max,N);    % angle vector [rad]
    else
        angles = 0;
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%   超声仿真   %%%%%%%%%%%%%%%%%%%%%%
    %% Source & sensor mask
    %% Based on the uff.probe we find the pixels in the domain that must work as source and sensors.
    trans_x = ((0:params.Nelements-1)-(params.Nelements-1)/2)*params.pitch; %以阵元中心位置为原点，每个阵元的位置
    % find the grid-points that match the element
    source_pixels={};
    element_sensor_index = {};
    n=1;  
    %% 阵元位置与仿真矩阵位置的匹配
    for m = 1:params.Nelements
        source_pixels{m} = find(abs(gridx+dx/2-trans_x(m))<=params.width/2 & gridz==0);
        element_sensor_index{m} = n:n+numel(source_pixels{m})-1;
        n=n+numel(source_pixels{m});
    end
    
    % sensor mask
    sensor.mask = zeros(size(medium.density));
    for m=1:params.Nelements
        sensor.mask(source_pixels{m}) = sensor.mask(source_pixels{m}) + 1;
    end
    
    % source mask
    source.u_mask=sensor.mask;
    
    % % % figure;
    % % % h = pcolor(gridx(1,:)*1e3,gridz(:,1)*1e3,source.u_mask); axis equal tight;
    % % % title('Source/Sensor mask')
    % % % set(h,'edgecolor','none');
    % % % set(gca,'YDir','reverse');
    % % % xlabel('x [mm]');
    % % % ylabel('z [mm]');
    
    %% Delay
    delay_values = zeros(params.Nelements, N);
    for n=1:N
        if angles(n) ~= 0
        delay_values(:,n) = -sqrt((trans_x-100*sin(angles(n))).^2+(zeros(size(trans_x))-100*cos(angles(n))).^2)./c0;
        delay_values(:,n) = delay_values(:,n)-min(delay_values(:,n));
        end
    end
    
    %% Calculation， launch the k-Wave calculation
    sensor_data = zeros(kgrid.Nt,params.Nelements*4,N); %% 改为*4
    disp('Launching kWave. This can take a while.');
    for n = 1:N
        delay = delay_values(:,n);
        denay = round(delay/kgrid.dt);
    %     seq(n).delay = seq(n).delay - cycles/f0/2;
        
        % offsets
        tone_burst_offset = [];
        for m=1:params.Nelements
            tone_burst_offset = [tone_burst_offset repmat(denay(m),1,numel(source_pixels{m}))];
        end
        source.ux = toneBurst(1/kgrid.dt, f0, cycles, 'SignalOffset', tone_burst_offset);   % create the tone burst signals
        source.uy = 0.*source.ux;
        source.u_mode ='additive'; 
        
        % set the input arguements: force the PML to be outside the computational
        % grid; switch off p0 smoothing within kspaceFirstOrder2D
        %input_args = {'DataCast', 'single','PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, 'Smooth', false};
        input_args = {'DataCast', 'gpuArray-single','PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, 'Smooth', false,'PlotScale','auto'};
        % run the simulation
        sensor_data(:,:,n) = gather(permute(kspaceFirstOrder2D(kgrid, medium, source, sensor, input_args{:}),[2 1]));
    end
    sensor_data(isnan(sensor_data)) = 0;
    
    %% Gather element signals, After calculaton we combine the signal recorded by the sensors according to the corresponding element
    element_data0 = zeros(numel(kgrid.t_array), params.Nelements, N);
    for m=1:params.Nelements
        if  ~isempty(element_sensor_index{m})
            element_data0(:,m,:)=bsxfun(@times,sqrt(1./kgrid.t_array).',trapz(kgrid.y(source_pixels{m}),sensor_data(:,element_sensor_index{m},:),2));
        end
    end
    element_data0(isnan(element_data0)) = 0;
    
    fs_simulation = 1 / kgrid.dt;      % 仿真信号采样率 [Hz]
    filtered_data = signal_pass_transducer_filter(element_data0, fs_simulation, params.fc, params.bandwidth, 'hanning');
    
    %% 重采样信号至 40 MHz
    element_data = resample(filtered_data, fs_simulation, params.fs);

    % % % plot(element_data(:,64,1));
    element_data(1:100,:) = 0;   % 除去激发的信号，如何更精准去除？
    element_data(isnan(element_data)) = 0;
    % % % plot((element_data(:,1:7:64)/max(element_data(:))+(1:10))',(0:size(element_data,1)-1)*kgrid.dt*1e6,'k')
    % % % set(gca,'XTick',1:10,'XTickLabel',int2str((1:7:64)'))
    % % % title('US RF signals')
    % % % xlabel('Element number'), ylabel('time (\mus)')
    % % % xlim([0 11])
    % % % axis ij square
    %%%  将探头收集到的超声RF数据及组织声阻抗率返回
    us_rf_data.rf_data = element_data;
    us_rf_data.acoustic_impedance = acoustic_impedance;
    
    %% Beamforming，To beamform we define a new (coarser) uff.linear_scan. We also define the processing pipeline and launch the beamformer
    % txdel = txdelay(0,100,param); % in s
    I = zeros(numel(z_imaging_axis), numel(x_imaging_axis));
    
    for n = 1:N
        txdel = delay_values;
        IQ = rf2iq(element_data, params.fs, params.fc);
        iq = IQ(:, params.Nelements);
        % % % plot(abs(iq))
        % % % set(gca,'YColor','none','box','off')
        % % % xlabel('time (\mus)')
        % % % title('I/Q signal')
        % % % % legend({'in-phase','quadrature'})
        % % % axis tight
        
        params.fnumber = [];
        [xi,zi] = meshgrid(x_imaging_axis, z_imaging_axis);
        
        bIQ = das(IQ, xi, zi, txdel, params);
        I(:,:) = bmode(abs(bIQ), 30); % log-compressed image, 这个值超声和光声应该不一样
    end
    %%% 将粗重建的超声图像返回
    us_img = I;
    
    % % % figure
    % % % imagesc(xi(1,:)*1e2, zi(:,1)*1e2, I)
    % % % colormap gray
    % % % title('US DAS image')
    % % % axis equal ij
    % % % set(gca,'XColor','none','box','off')
    % % % c = colorbar;
    % % % c.YTick = [0 255];
    % % % c.YTickLabel = {'-40 dB','0 dB'};
    % % % ylabel('[cm]')
    
    
    %%%%%%%%%%%%%%%%%%%%%%%   光声仿真   %%%%%%%%%%%%%%%%%%%%%%
    %% 介质与超声仿真时一样，不再重新定义
    %% define the time vector depending on the CFL number, the size of the domain and the mean speed of sound.
    cfl = 0.3;
    t_end = sqrt(grid_depth.^2+grid_depth.^2)/mean(medium.sound_speed(:));  % 矩阵对角线单程
    kgrid.makeTime(medium.sound_speed, cfl, t_end);
    
    % source mask，将加载的图像作为声源，初始声压
    source_PA.p_mask = single(logical(p0-1));
    source_PA.p0 = optical_absorption * 1.*abs((p0-1));  %% 把不同组织对光吸收差异（光声压差异）考虑进去
    % % % figure;
    % % % h = pcolor(gridx(1,:)*1e3,gridz(:,1)*1e3,source_PA.p0); axis equal tight;
    % % % title('Source/Sensor mask')
    % % % set(h,'edgecolor','none');
    % % % set(gca,'YDir','reverse');
    % % % xlabel('x [mm]');
    % % % ylabel('z [mm]');
    
    %% Calculation， launch the k-Wave calculation
    sensor_data_PA = zeros(kgrid.Nt,params.Nelements*4,N);
    disp('Launching kWave. This can take a while.');
    % set the input arguements: force the PML to be outside the computational
    % grid; switch off p0 smoothing within kspaceFirstOrder2D
    %input_args = {'DataCast', 'single', 'PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, 'Smooth', false};
    input_args = {'DataCast', 'gpuArray-single','PMLInside', false, 'PMLSize', PML_size, 'PlotPML', false, 'Smooth', false,'PlotScale','auto'};
        % run the simulation
    sensor_data_PA(:,:,n) = gather(permute(kspaceFirstOrder2D(kgrid, medium, source_PA, sensor, input_args{:}),[2 1]));
    sensor_data_PA(isnan(sensor_data_PA)) = 0;
    
    %% Gather element signals，After calculaton we combine the signal recorded by the sensors according to the corresponding element
    element_data_PA0 = zeros(numel(kgrid.t_array),params.Nelements,N);
    for m = 1:params.Nelements
        if  ~isempty(element_sensor_index{m})
            element_data_PA0(:,m,:)=bsxfun(@times,sqrt(1./kgrid.t_array).',trapz(kgrid.y(source_pixels{m}),sensor_data_PA(:,element_sensor_index{m},:),2));
        end
    end
    element_data_PA0(isnan(element_data_PA0)) = 0;
    
    %% 对采集的光声信号进行滤波，滤波参数与超声时一样
    filtered_data_PA = signal_pass_transducer_filter(element_data_PA0, fs_simulation, params.fc, params.bandwidth, 'hanning');
    
    %% 重采样信号至 40 MHz
    element_data_PA = resample(filtered_data_PA, fs_simulation, params.fs);
    element_data_PA(isnan(element_data_PA)) = 0;
    %%%  将探头收集到的光声RF数据返回
    pa_rf_data.rf_data = element_data_PA;
    pa_rf_data.optical_absorption = optical_absorption;
    % % % plot((element_data_PA(:,1:7:64)/max(element_data_PA(:))+(1:10))', (0:size(element_data_PA,1)-1)*kgrid.dt*1e6,'k')
    % % % set(gca,'XTick',1:10,'XTickLabel',int2str((1:7:64)'))
    % % % title('PA RF signals')
    % % % xlabel('Element number'), ylabel('time (\mus)')
    % % % xlim([0 11])
    % % % axis ij square
    
    %% Beamforming，To beamform we define a new (coarser) uff.linear_scan. We also define the processing pipeline and launch the beamformer
    I_PA = zeros(numel(z_imaging_axis), numel(x_imaging_axis));
    
    txdel = delay_values;
    IQ_PA = rf2iq(element_data_PA, params.fs, params.fc);
    iq_PA = IQ_PA(:, params.Nelements);
    % % % plot(abs(iq_PA))
    % % % set(gca,'YColor','none','box','off')
    % % % xlabel('time (\mus)')
    % % % title('I/Q signal')
    % % % % legend({'in-phase','quadrature'})
    % % % axis tight
    
    params.fnumber = [];
    [xi,zi] = meshgrid(x_imaging_axis, z_imaging_axis);
    params_PA = params;
    params_PA.passive = true;
    bIQ_PA = das(IQ_PA, xi, zi, txdel, params_PA);
    
    I_PA(:,:) = bmode(abs(bIQ_PA), 20); % log-compressed image
    
    %%%  将 粗重建的光声图像存盘以备后续网络使用
    pa_img = I_PA;
    
    % % % figure
    % % % imagesc(xi(1,:)*1e2,zi(:,1)*1e2,I_PA)
    % % % colormap gray
    % % % title('PA DAS image')
    % % % axis equal ij
    % % % set(gca,'XColor','none','box','off')
    % % % c = colorbar;
    % % % c.YTick = [0 255];
    % % % c.YTickLabel = {'-40 dB','0 dB'};
    % % % ylabel('[cm]')
end

%%% 信号通过换能器的滤波器函数, TODO 输入信号也可能是一维的
%%
%% input_signal: 输入信号
%% sampling_frequency: 信号采样频率
%% frequency_center: 换能器中心频率
%% band_width: 换能器带宽，百分比形式[%]
%% filter_window: 滤波器窗口类型
function output_signal = signal_pass_transducer_filter(input_signal, sampling_frequency, frequency_center, band_width, filter_window)
    band_width = band_width / 100.0;   % 带宽百分比
    % 计算带通滤波器的截止频率
    f_low = frequency_center * (1 - band_width / 2);
    f_high = frequency_center * (1 + band_width / 2);
    
    % 设计带通滤波器（FIR 滤波器）
    filter_order = 100; % 滤波器阶数，可根据需要调整
    b = fir1(filter_order, [f_low, f_high] / (sampling_frequency / 2), 'bandpass');
    
    % 对信号进行滤波
    filtered_data = zeros(size(input_signal)); % 预分配内存
    for i = 1:size(input_signal, 2)
        filtered_data(:, i) = filtfilt(b, 1, input_signal(:, i)); % 双向滤波，保持相位
    end 

    output_signal = filtered_data;
end

%%% 重采样信号至目标采样率
%%
%% input_signal:  输入信号
%% sampling_rate: 输入信号的采样频率
%% target_sampling_rate: 目标采样频率
function output_signal = resample(input_signal, sampling_rate, target_sampling_rate) 
    % 获取原始采样时间轴
    t_original = (0:size(input_signal, 1) - 1) / sampling_rate;
    
    % 生成目标采样时间轴
    t_target = (0:1/target_sampling_rate:(t_original(end))); 
    
    % 初始化存储重采样信号的矩阵
    resampled_data = zeros(length(t_target), size(input_signal, 2));
    
    % 对每列信号进行重采样
    for i = 1:size(input_signal, 2)
        resampled_data(:, i) = interp1(t_original, input_signal(:, i), t_target, 'linear');
    end
    output_signal = resampled_data;
end
    